********************************************************************************
* LPDENSITY STATA PACKAGE -- lpdensity
* Authors: Matias D. Cattaneo, Michael Jansson, Xinwei Ma
********************************************************************************
*!version 1.0 14-Jul-2019 

capture program drop lpdensity

program define lpdensity, eclass
syntax 	varlist(max=1) 			///
								///
		[if] [in] [, 			///
								///
		GRId(varname) 			///
		BW(string) 			///
		P(integer 2) 			///
		Q(integer -1) 			///
		V(integer 1) 			///
		BWSelect(string) 		///
		KERnel(string) 			///
		CWeights(varname) 		///
		PWeights(varname) 		///
		SCAle(real 1) 			///
		GENvars(string) 		///
		Level(real 95) 			///
		SEParator(integer 5)	///
								///
		PLot					///
		GRAph_options(string)	///
		]

marksample touse

********************************************************************************
**** error check: main variable
local x "`varlist'"
cap confirm numeric variable `x'
if _rc {
	di as err `"`grid' is not a numeric variable"'
	exit 198
}
else {
	qui count if `x' != . & `touse'
	local n = r(N)
	if (`n' == 0) {
		di as err `"`x' has length 0"'
		exit 198
	}
}

********************************************************************************
**** error check: grid()
if ("`grid'" == "") {
	tempvar temp_grid	
	qui pctile `temp_grid' = `x' if `touse', nq(20)
	local ng = 19
} 
else {
	cap confirm numeric variable `grid'
	if _rc {
		di as err `"grid(): `grid' is not a numeric variable"'
		exit 198
	}
	else {
		if ("`x'" == "`grid'") {
			qui count if `grid' != . & `touse'
		}
		else {
			qui count if `grid' != .
		}
		
		local ng = r(N)
		//disp `ng'
		if (`ng' == 0) {
			di as err `"grid(): `grid' has length 0"'
			exit 198
		}
	}
	local temp_grid "`grid'"
}

********************************************************************************
**** error check: p() q() v()
if (`p' < 0 | `p' > 7) {
	di as err `"p(): has to be integer between 0 and 7"'
	exit 198
}
if (`q' == -1) {
	local q = `p' + 1
}
if (`q' < `p' | `q' > 7) {
	di as err `"q(): has to be integer between `p' and 7"'
	exit 198
}
if (`v' < 0 | `v' > `p') {
	di as err `"v(): has to be integer between 0 and `p'"'
	exit 198
}

********************************************************************************
**** error check: bwselect()
if ("`bwselect'" == "") {
	local bwselect = "mse-dpi"
} 
else {
	if ("`bwselect'" != "mse-dpi" & "`bwselect'" != "imse-dpi" & "`bwselect'" != "mse-rot" & "`bwselect'" != "imse-rot") {
		di as err `"bwselect(): incorrectly specified: options(mse-dpi, imse-dpi, mse-rot, imse-rot)"'
		exit 198
	}
}

********************************************************************************
**** error check: kernel()
if ("`kernel'" == "") {
	local kernel = "triangular"
} 
else {
	if ("`kernel'" != "triangular" & "`kernel'" != "uniform" & "`kernel'" != "epanechnikov") {
		di as err `"kernel(): incorrectly specified: options(triangular, uniform, epanechnikov)"'
		exit 198
	}
}

********************************************************************************
**** error check: scale()
if (`scale' < 0) {
	di as err `"scale(): incorrectly specified: should be between 0 and 1"'
	exit 198
}

********************************************************************************
**** error check: level()
if (`level' <= 0 | `level' >= 100) {
	di as err `"level(): incorrectly specified: should be between 0 and 100"'
	exit 198
}

********************************************************************************
**** error check: cweights()
if ("`cweights'" == "") {
	tempvar temp_cweights
	qui gen `temp_cweights' = `x' * 0 + 1
}
else {
	cap confirm numeric variable `cweights'
	if _rc {
		di as err `"cweights(): `cweights' is not a numeric variable"'
		exit 198
	}
	else {
		qui count if `cweights' != . & `touse'
		local nCw = r(N)
		if (`nCw' != `n') {
			di as err `"cweights(): `cweights' has different length as `x'"'
			exit 198
		}
	}
	local temp_cweights "`cweights'"
}

********************************************************************************
**** error check: pweights()
if ("`pweights'" == "") {
	tempvar temp_pweights
	qui gen `temp_pweights' = `x' * 0 + 1
}
else {
	cap confirm numeric variable `pweights'
	if _rc {
		di as err `"pweights(): `pweights' is not a numeric variable"'
		exit 198
	}
	else {
		qui count if `pweights' != . & `touse'
		local nPw = r(N)
		if (`nPw' != `n') {
			di as err `"pweights(): `pweights' has different length as `x'"'
			exit 198
		}
	}
	local temp_pweights "`pweights'"
}
 
********************************************************************************
**** error check: bw()
if ("`bw'" == "") {
	tempvar temp_bw
	tempvar temp_lpbwdensity
	qui lpbwdensity `x' if `touse', grid(`temp_grid') p(`p') v(`v') bwselect(`bwselect') kernel(`kernel') cweights(`temp_cweights') pweights(`temp_pweights') genvars(`temp_lpbwdensity') separator(1)
	qui gen `temp_bw' = `temp_lpbwdensity'_bw
	qui capture drop `temp_lpbwdensity'_grid `temp_lpbwdensity'_bw `temp_lpbwdensity'_nh
}
else {

	cap confirm numeric variable `bw'
	if _rc {
		// check if variable exists
		capture confirm variable `bw'
		if (!_rc) {
			di as err `"bw(): `bw' is not a numeric variable"'
			exit 198
		} 
		else {
			tempvar temp_bw			
			qui gen `temp_bw' = "`bw'" if `temp_grid' != .
			capture destring `temp_bw', replace
			cap confirm numeric variable `temp_bw'
			if _rc {
				di as err `"bw(): `bw' is not numeric"'
				exit 198
			}
		}
	}
	else {
		if ("`x'" == "`grid'") {
			qui count if `bw' != . & `touse'
		}
		else {
			qui count if `bw' != .
		}
		
		local nb = r(N)
		if (`nb' != `ng') {
			di as err `"bw(): `bw' has different length as grid()"'
			exit 198
		}
		local temp_bw "`bw'"
		local bwselect = "user provided"
	}
}

********************************************************************************
**** error check: separator()
if (`separator' <= 1) {
	local separator = 1
}

********************************************************************************
**** temporaty varaibles for plotting
if ("`plot'" != "" & "`genvars'" == "") {
	tempvar plot_grid
	tempvar plot_f
	tempvar plot_cil
	tempvar plot_cir
}

********************************************************************************
**** MATA
tempvar temp_touse
qui gen `temp_touse' = `touse'
mata{
	x = st_data(., "`x'", "`temp_touse'") //; x
	
	if ("`x'" == "`grid'") {
		grid     = st_data(., "`temp_grid'"	, "`temp_touse'") //; grid
		bw       = st_data(., "`temp_bw'"	, "`temp_touse'") //; bw
	}
	else {	
		grid     = st_data(., "`temp_grid'"	, 0) //; grid
		bw       = st_data(., "`temp_bw'"	, 0) //; bw
	}

	cweights = st_data(., "`temp_cweights'"	, "`temp_touse'") //; cweights
	pweights = st_data(., "`temp_pweights'"	, "`temp_touse'") //; pweights
	
	Result = lpdensity_fn(x, grid, bw, `p', `q', `v', "`kernel'", cweights, pweights, `scale', 1-`level'/100)
	
	genvars = st_local("genvars")
	if (genvars != "") {
		(void) 	st_addvar("double", 	invtokens((genvars, "grid"), 		"_"))
				st_store(Result[., 10], invtokens((genvars, "grid"), 		"_"), Result[., 1])
					
		(void) 	st_addvar("double", 	invtokens((genvars, "bw"), 			"_"))
				st_store(Result[., 10], invtokens((genvars, "bw"), 			"_"), Result[., 2])
					
		(void) 	st_addvar("double", 	invtokens((genvars, "nh"), 			"_"))
				st_store(Result[., 10], invtokens((genvars, "nh"), 			"_"), Result[., 3])
					
		(void) 	st_addvar("double", 	invtokens((genvars, "f_p"), 		"_"))
				st_store(Result[., 10], invtokens((genvars, "f_p"), 		"_"), Result[., 4])
					
		(void) 	st_addvar("double", 	invtokens((genvars, "se_p"), 		"_"))
				st_store(Result[., 10], invtokens((genvars, "se_p"), 		"_"), Result[., 6])
					
		(void) 	st_addvar("double", 	invtokens((genvars, "CI_l"), 		"_"))
				st_store(Result[., 10], invtokens((genvars, "CI_l"), 		"_"), Result[., 8])
					
		(void) 	st_addvar("double", 	invtokens((genvars, "CI_r"), 		"_"))
				st_store(Result[., 10], invtokens((genvars, "CI_r"), 		"_"), Result[., 9])
					
		if (`q' > `p') {
		(void) st_addvar("double", 		invtokens((genvars, "f_q"), 		"_"))
				st_store(Result[., 10], invtokens((genvars, "f_q"), 		"_"), Result[., 5])
		
		(void) st_addvar("double", 		invtokens((genvars, "se_q"), 		"_"))
				st_store(Result[., 10], invtokens((genvars, "se_q"), 		"_"), Result[., 7])
		}
	}
	else if ("`plot'" != "") {
		(void) 	st_addvar("double", 	"`plot_grid'")
				st_store(Result[., 10], "`plot_grid'", Result[., 1])
				
		(void) st_addvar("double", 		"`plot_f'")
				st_store(Result[., 10], "`plot_f'", Result[., 4])
				
		(void) st_addvar("double", 		"`plot_cil'")
				st_store(Result[., 10], "`plot_cil'", Result[., 8])
				
		(void) st_addvar("double", 		"`plot_cir'")
				st_store(Result[., 10], "`plot_cir'", Result[., 9])
	}
	st_matrix("Result", Result[., 1..9])
}

********************************************************************************
**** generate variable labels
if ("`genvars'" != "") {
    label variable `genvars'_grid 	"lpdensity: grid point"
	label variable `genvars'_bw 	"lpdensity: bandwidth"
	label variable `genvars'_nh 	"lpdensity: effective sample size"
	label variable `genvars'_f_p 	"lpdensity: point estimate with pol. order `p'"
	label variable `genvars'_se_p 	"lpdensity: standard error for f_p"
	label variable `genvars'_CI_l 	"lpdensity: left level-`level' CI, conventional"
	label variable `genvars'_CI_r 	"lpdensity: right level-`level' CI, conventional"
	if (`q' > `p') {
	label variable `genvars'_f_q 	"lpdensity: point estimate with pol. order `q'"
	label variable `genvars'_se_q 	"lpdensity: standard error for f_q"
	label variable `genvars'_CI_l 	"lpdensity: left level-`level' CI, robust bias-corrected"
	label variable `genvars'_CI_r 	"lpdensity: right level-`level' CI, robust bias-corrected"
	}	
}

********************************************************************************
**** display
disp ""
disp "Local Polynomial Density Estimation and Inference." 
disp ""

disp in smcl in gr "{lalign 1: Sample size                              (n=)    }" _col(19) in ye %15.0f `n'
disp in smcl in gr "{lalign 1: Polynomial order for point estimation    (p=)    }" _col(19) in ye %15.0f `p'
if (`v' == 1) {
disp in smcl in gr "{lalign 1: Density function estimated               (v=)    }" _col(19) in ye %15.0f `v'
}
else if (`v' == 0) {
disp in smcl in gr "{lalign 1: Distribution function estimated          (v=)    }" _col(19) in ye %15.0f `v'
} 
else {
disp in smcl in gr "{lalign 1: Order of derivative estimated            (v=)    }" _col(19) in ye %15.0f `v'
}
disp in smcl in gr "{lalign 1: Polynomial order for confidence interval (q=)    }" _col(19) in ye %15.0f `q'
disp in smcl in gr "{lalign 1: Kernel function                                  }" _col(19) in ye "{ralign 15: `kernel'}"
if (`scale' < 1) {
disp in smcl in gr "{lalign 1: Scaling factor                                   }" _col(19) in ye %15.2f `scale'
}
disp in smcl in gr "{lalign 1: Bandwidth selection method                       }" _col(19) in ye "{ralign 15: `bwselect'}"
disp ""

disp in smcl in gr "{hline 72}"
if (`q' > `p') {
disp in smcl in gr "{ralign 4: }" _col(4) "{ralign 10: }" _col(14) "{ralign 10: }" _col(24) "{ralign 8: }" _col(32) "{ralign 10: Point}" _col(42) "{ralign 10: Std.}" _col(52) "{ralign 20: Robust B.C.}" _col(72)
} 
else {
disp in smcl in gr "{ralign 4: }" _col(4) "{ralign 10: }" _col(14) "{ralign 10: }" _col(24) "{ralign 8: }" _col(32) "{ralign 10: Point}" _col(42) "{ralign 10: Std.}" _col(52) "{ralign 20: Conventional}" _col(72)
}
disp in smcl in gr "{ralign 4: }" _col(4) "{ralign 10: grid}" _col(14)  "{ralign 10: bw}" _col(24) "{ralign 8: Eff.n}" _col(32) "{ralign 10: Est.}" _col(42) "{ralign 10: Error}" _col(52) "{ralign 20: `level'% Conf. Interval}" _col(72)
disp in smcl in gr "{hline 72}"

forvalues i = 1(1)`ng' {
disp in smcl in gr %4.0f `i' _col(4) in ye %10.4f Result[`i', 1] _col(14)  %10.4f Result[`i', 2] _col(24) %8.0f Result[`i', 3] _col(32) %10.4f Result[`i', 4] _col(42) %10.4f Result[`i', 6] _col(52) ///
    in ye %10.4f Result[`i', 8] _col(62) in ye %10.4f Result[`i', 9] _col(72)
	if (`i' != `ng' & `separator' > 1 & mod(`i', `separator') == 0) {
		disp in smcl in gr "{hline 72}"
	}
}
disp in smcl in gr "{hline 72}"

********************************************************************************
**** plot
if ("`plot'" != "" & "`genvars'" != "") {
	local plot_grid	"`genvars'_grid"
	local plot_f	"`genvars'_f_p"
	local plot_cil	"`genvars'_CI_l"
	local plot_cir 	"`genvars'_CI_r"
}

if ("`plot'" != "") {
	if (`"`graph_options'"'=="" ) local graph_options = `"title("lpdensity (p=`p', q=`q', v=`v')", color(gs0)) xtitle("`x'") ytitle("")"'
	twoway 	(rarea `plot_cil' `plot_cir' `plot_grid', sort color(gs11)) ///
			(line `plot_f' `plot_grid',  lcolor(black) sort lwidth(medthin) lpattern(solid)),  ///
			legend(cols(2) order(2 "point estimate" 1 "`level'% C.I." )) `graph_options'
}

********************************************************************************
**** ereturn

ereturn clear
ereturn scalar N = `n'
ereturn scalar p = `p'
ereturn scalar q = `q'
ereturn scalar v = `v'
ereturn local kernel = "`kernel'"
ereturn local bwselect = "`bwselect'"
ereturn scalar scale = `scale'
ereturn scalar level = `level'
matrix colnames Result = grid bw nh f_p f_q se_p se_q CI_l CI_r
ereturn matrix result = Result

********************************************************************************
**** end

mata: mata clear

end
	
